Dark matter clustering: a simple renormalization group approach 



Patrick McDonald^ 

Canadian Institute for Theoretical Astrophysics, University of Toronto, Toronto, ON M5S 3H8, Canada 

(Dated: February 5, 2008) 

I compute a renormalization group (RG) improvement to the standard beyond-linear-order Eule- 
rian perturbation theory (PT) calculation of the power spectrum of large-scale density fluctuations 
in the Universe. At z = 0, for a power spectrum matching current observations, lowest order RGPT 
appears to be as accurate as one can test using existing numerical simulation-calibrated fitting 
formulas out to at least k ~ 0.3 h Mpc -1 ; although inaccuracy is guaranteed at some level by ap- 
proximations in the calculation (which can be improved in the future). In contrast, standard PT 
breaks down virtually as soon as beyond-linear corrections become non-negligible, on scales even 
larger than k = 0.1 /iMpc -1 . This extension in range of validity could substantially enhance the 
, usefulness of PT for interpreting baryonic acoustic oscillation surveys aimed at probing dark energy, 

■ for example. I show that the predicted power spectrum converges at high k to a power law with 

index given by the fixed-point solution of the RG equation. I discuss many possible future directions 
£\J , for this line of work. The basic calculation of this paper should be easily understandable without 

any prior knowledge of RG methods, while a rich background of mathematical physics literature 
exists for the interested reader. 
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Statistics of the large-scale density fluctuations in the Universe provide an invaluable source of information on the 
very early Universe, where the initial density perturbations that grow by gravitational instability were created, and 
on the material content and gravitational laws of the Universe, which determine the growth of structure and how it 
looks to us from a distance [l[ . The standard practical toolbox we use for interpreting observations consists of linear 
theory on very large scales or at very early times [2(], and everywhere else numerical simulations which compute the 
fully non-linear evolution of a realization of the density field in a patch of Universe. Fitting formulas, motivated by 
physical intuition but ultimately heuristic, are often used to more efficiently interpolate between simulation results 



Beyond linear order perturbation theory (PT) has been considered for a long time [H, foL ItI l9l, fiol. ITlj ■ but does not 



currently play an important role in our interpretation of precision observations (see jl2j for a thorough review of large- 
scale structure perturbation theory). The desire to measure dark energy properties using observations of baryonic 



scale structure perturbation theory), ine desire to measure dark energy properties using observations ol baryonic 
acoustic oscillation features [H EI EE ES E3, El El, HI HH H3, [H, HI HE IM H3 , which fal1 more or less exactly 
in the range of scales where higher order perturbation theory could be most useful, gives this line of work much more 
pressing practical relevance than it has had in the past [281 ]. More traditional uses of galaxy clustering to measure 
; I ■ cosmological parameters, long reliant on the assumption that galaxy power is simply proportional to linear mass 
power, have also reached the point where weakly non-linear effects are a substantial limitation [29L [3fjj. Additionally, 
other probes of large-scale structure, like the Lya forest [3l|, l32|, weak lensing [33| . galaxy cluster/Sunyaev-Zel'dovich 
effect (SZ) measurements [34 |. and possibly future 21cm surveys [35| . could all benefit from improved computational 
techniques. 

The problem with higher order perturbation theory is that it does not work very well in exactly the regime where 
it could currently be most useful: at low redshift and wavenumber k ~ 0.1 /iMpc" 1 , where galaxy clustering is 
significantly but still only weakly non-linear [36| . Recently, [37l . |38| proposed a renormalization technique aimed at 
improving PT calculations. They re-sum an infinite series of perturbation theory terms (Feynman diagrams) which 
determine the memory of perturbations to initial conditions as a function of scale (see also [39]). This method is 
apparently computationally difficult because they have not yet computed a power spectrum using it. The renormal- 
ization group (RG) method I present here is different, but surely related. There has been some work toward using 
RG techniques to model large-scale structure [13, EH, S3] ■ The methods used are different than mine, and have not 
yet seen much practical application for comparison with observations. 
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My implementation of the RG is inspired more by its use for dealing with secular divergences in perturbative 
solutions of differential equations [H, [44[ than by the probably better known applications in high energy and statistical 
physics [H, H(| IU EH , although all of these uses are related. The basic idea is to identify a quantity in the problem 
that is not directly observable and modify it to remove breakdowns in the perturbation theory, leading to a well- 
behaved renormalized observable. In this paper I renormalize the initial (linear theory) power spectrum to remove 
the relative (not actually infinite for realistic power spectra) divergence of the lowest order correction to linear theory. 



More detailed mathematical clarification and explanation of the meaning of the method of 1431. 14411 . in terms of the 
theory of envelopes, can be found in [49|,|50|. Even more detailed treatments can be found in IHJII. These papers 



should be useful to the reader who is unsatisfied with my presentation, or wants to extend it. During the last several 
years, [53|, H3, [55|, H^, H3, [H, H^] have extended and generalized this dynamical renormalization group method to 
compute the real time evolution of expectation values of quantum fields. 

The rest of the paper is as follows: In |Jl]l derive a simple formula for the RG- improved power spectrum of mass 
density fluctuations. In ^IIII I show numerical results and compare them to standard PT and fitting formulas calibrated 
using simulations. Finally, in §IVII discuss possible future directions for this kind of work. 

II. CALCULATION 

In this section I present the basic calculation. In All Al I discuss the time evolution equations for cold dark matter, 
which are the starting point for perturbation theory. In EIIIBI I review standard Eulerian perturbation theory. Finally, 
in mi CI I apply the renormalization group improvement. 

A. Evolution equations 

We are interested in the statistics of the mass density field, <5(x,t) = p(x,r)/p — 1, with Fourier transform S(k, r) = 
J g? 3 x exp(ik • x) <5(x, r), where x is the comoving position and r = J dt/a is the conformal time, with a = 1/(1 + z) 
the expansion factor. In particular, I will compute the power spectrum, P(k,r), defined by 

(S(k, r)<5(k',r)) = (2ir) 3 5 D (k + k')P(k, r) (1) 

For now I will assume only cold dark matter, i.e., collisionlcss particles with insignificant initial velocity dispersion. 
The exact evolution is is described by the Vlasov equation [j| 

df 1 

J + p.V/- fl mVfV p / = 0, (2) 

or am 

with 

V 2 </> = 4 7T G a 2 p S , (3) 

where /(x, p, r) is the particle density at phase-space position (x, p), m is the particle mass (which plays no role in 
the final results), and p = a mv (v here is a particle's peculiar velocity, not to be confused with the mean peculiar 
velocity of all particles at some point in space, used everywhere else). Except when otherwise indicated V = d/dx. 
The density field is obtained by averaging the distribution function over momentum: 

p(x, t) = ma- 3 J d 3 p f(x, p, r) , (4) 

and the bulk (mean) velocity and higher moments of the velocity distribution e.g., the dispersion of particle velocities 
around their bulk velocity, can be similarly obtained by multiplying the distribution function by any number of p's 
(e.g., one to obtain bulk velocity) before integrating over p. As discussed by [5[, taking moments of the Vlasov 
equation with respect to momentum leads to a hierarchy of evolution equations for these quantities. Dropping the 
higher moments of the velocity distribution gives the usual hydrodynamic equations for density and bulk velocity: 

g + V .[(l + 5)v]=0, (5) 

and 

dv 

F +Hv|(vV)v = -V^, (6) 



where TL = dlna/dr = aH with H the usual Hubble parameter 
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B. Standard Eulerian perturbation theory 



Perturbation theory consists of writing the density and velocity fields as a series of terms of at least formally 
increasing order of smallness, i.e., <5 = Si + #2 + #3 + •••• The evolution equations are solved order- by-order, with lower 
order solutions appearing as sources in the higher order equations so that S n is of order <5™ [l2| . The power spectrum 
for Gaussian initial conditions is given by 

(*(k)a(k')> = (fcOWk')) + (SxQLfaW) + (* 3 (k)* 1 (k')> + (<5 2 (k)<5 2 (k')) + ... (7) 

where no terms 3rd order in 5% appear because the expectation value of any term cubic in a Gaussian field is zero. 
At this point I assume an Einstein-de Sitter Universe for simplicity, and define 



P(k,r) = D 2 (r)P n (k) + D\t) [P 13 (k) + P 22 (k)] + 



(8) 



where D(t) = Si (r) /initial is the linear theory growth factor. The Einstein-de Sitter assumption is needed to 
avoid more complicated time dependence of Pi 3 and P22, but the real Universe is of course not Einstein-de Sitter. 
Fortunately, this is not a significant problem because, to percent level accuracy [28|, HO, [Ml H3, [fill H3, [65|, [6||, the 
effect of changing the background model can be included by simply using the correct linear growth factor in Eq. (JSJ). I 
will sometimes refer to Pi 3 and P22 as 2nd order, meaning in the initial power spectrum amplitude, not to be confused 
with the fact that they are 4th order in Si and require calculating the evolution of S to 3rd order, i.e., S3. [67j derived 
the following useful form of the equations for Pi3(fc) and Pa2(A0 : 



fc 3 Pn(fc) 
252 (2tt) 2 Jo 



drPn(kr) 
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158 + 100r 2 - 42r 4 



and 



P22W 



k 3 



98 (2tt) z Jo 



drPn (kr) / dxPu k (l 



2rx 



1/2 



1) (7r 2 + 2)ln 



(3r + 7x- Wrx 2 ) 
(1 + r 2 - 2ra;) 2 



1 



(9) 



(10) 



Note that, due to many cancellations, these terms are not as divergent as they might appear at first glance, their sum 
being convergent at high k for power law Pn(fc) with n = d\nP/d\nk < —1 and at low k when n > —3 [67j |. Their 
sum is zero for n ~ —1.4 [l2j], a fact that will have interesting consequences for our RG calculation. 



C. Renormalization group improvement 

The problem with the standard calculation outlined in ^IIBl which leads us to renormalization, is that the 2nd 
term in Eq. ([5]) diverges relative to the first (although it does not literally become infinite for realistic power spectra), 
at increasingly large scales (small k) as time progresses. This divergence is not exactly unphysical - we expect non- 
linearities to become important during gravitational collapse - however, the basic premise of perturbation theory is 
violated when higher order terms become large. There is no a priori reason to expect the results to be accurate. I 
employ a renormalization group calculation to cure this divergence. I start in All C ll with a quick derivation motivated 
by the method of [43j . Then I discuss the mathematical background, physical interpretation, approximations, and 
limitations of the method in §11 C 21 



1. Simple calculation 

To simplify the presentation of the calculation, I rewrite Eq. ((8|) in a more compact, schematic form, 

P(k,r)^P L (k)+A [Pl](k) , (11) 

where A cx D 2 (t), P(fc, t) = P(k,r)/A, Ph{k) is the initial condition power, and [P|](/c) = Pis(k) + p22(&) is the 
higher order correction term, which is quadratic in Pt(fc), as given by Eqs. (0) and (|10|) . I now follow the method 
described in (43| to deal with the problem of the 2nd order term becoming large (deferring a detailed explanation of 
the meaning of this method to mi C 2[) . I rewrite A = A — A+ + A*, where A+ is just an arbitrary constant, so that 



P(k, t) = P L (k) + A* [Pi] (k) + (A - A,) [Pi] (k) 



(12) 
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I then absorb the constant A+ [P|] (k) into the initial conditions, Pl(£;), defining the renormalized initial conditions 
P*(fc, j4*), so that 

P(k) = P*(k, A) + (A - A) [P 2 ] (fc) , (13) 

where I can freely replace P^ with P* in the 2nd term because the change is formally 3rd order (in the initial power 
spectrum amplitude). The observable power spectrum P(k) obviously must be independent of the completely artificial 
constant A*, so I impose dP(k) /dA± = and find 

dP f A ^ ] - [PR (*) , (14) 

where I have dropped the derivative of [P 2 ] (k) because it is higher order (this is seen to be self-consistent from Eq.[T4l 
itself). This equation is easy enough to solve numerically, given P*(fc, A*) at some value of A*. Since A+ is completely 
arbitrary, I will in the end choose A+ = A, removing the 2nd term in Eq. (fT"3]) entirely. The result for P(fc,r) is then 
simply P*(k,A± = A), obtained by solving Eq. (fl"4"|) . To reproduce the linear theory result at very early times, the 
initial conditions for the solution of Eq. JTJJ must be P*(fc, A* = 0) — Pi(fc). 



2. Explanation 

The method of [43| and related papers was, b y th e authors own admission, somewhat lacking in rigorous motivation. 
This situation was clarified by others (49l . [Hoi. l5ll. [52J; however, their explanations may be opaque. A very simple 
example should serve to explain and clarify the method. Consider the differential equation 

S = S + S 2 . (15) 

The exact solution is easy to obtain, but suppose one did not know that and solved Eq. (j!5[) using perturbation 
theory for small initial 8, as in the cosmological situation. The solution to the linearized equation is Si — <?ie*. The 
equation for the first non-linear correction, S2, is then 

52=S 2 +gje 2t (16) 

with solution 

62 =g2e t +gfe 2t . (17) 

The approach generally followed in the cosmological case is to set the term that is the solution to the homogeneous 
equation, c^e*, to zero, i.e., to assume that at asymptotically early times 5 — ► <5i; however, this is a choice, not 
necessary. The key observation here is that solving two differential equations has produced two parameters in the 
solution, g\ and g 2 , but we only need one to satisfy the physical boundary conditions. The RG method exploits and 
fixes the otherwise arbitrary extra parameter of the perturbative solution (note that this approach gives exactly the 
same result as the approach of "splitting" the linear solution into two pieces, one of which is considered to be second 
order and used to cancel the non-linear correction). 

For the standard choice g 2 = 0, 62 inevitably becomes larger than Si, invalidating the perturbation theory. Note, 
however, that using g2 we can actually guarantee that at any one chosen time, £*, S% (£*) = 0, and very near t+, 
82 -C 6\. This requires g2 = — sfe 4 *, which leads to full solution 

S^gie'+gje'ie'-e^). (18) 

Clearly the value of gi must depend on the choice of i.e., gi = <7i (£*.), for a fixed physical situation (fixed full 
solution (5(f)); however, nothing would be accomplished by simply fixing gi (f*) to match the initial conditions, because 
this generally requires using the perturbative solution outside its range of validity. The RG method is to determine 
<7i(f*) by enforcing the idea that the physical solution should be independent of the arbitrary time f*, while using the 
perturbative solution only for f infinitcsimally close to f*, where it should be valid, i.e., 



dS(t) 



dt* 

or 



= 0=^e*-«^ (19) 
u=t at 



(20) 



5 



The exact solution to this equation is <?i(t) = c/(l — ce*) where c is a constant which will allow matching of the 
physical initial conditions. Finally, referring back to Eq. (fT5|) , evaluated at = i, gives the final solution 

*(*) = T~~t- (21) 
1 — ce 1 

While such good fortune would generally not be expected, it happens that this is the exact solution to the original 
differential equation, Eq. (fT5)) . 

The intuitive picture one should have of this calculation is of stepping forward in time slightly using the perturbative 
solution centered at the present time, then taking the result and using it as the initial condition for a perturbative 
solution around the new time, followed by another step, and so on. [49j showed that this method gives the envelope 
equation for the set of solutions obtained by applying the perturbation theory near various times t* (the envelope is 
the curve which is tangent to each of the different solutions). 

To make the analogy with original case of the power spectrum clear, note that changing variables to A = exp (t), 
A+ = exp(t*), and 6 = 8/ A in Eq. (fig)) produces an equation with the same form as Eq. (Tj"3")) . The rest of the 
derivation follows identically until, unfortunately, Eq. (TTJJ is not analytically solvable because P(k) represents an 
infinite number of variables rather than one. 

While the resulting equations look similar, the skeptical reader may still question the relation between the derivation 
of Eq. (fT4|) for the power spectrum and Eq. (j20|) for a single variable governed by a differential equation, e.g., S would be 
most naturally be identified with the cosmological density field, not the power spectrum. To complete the connection 
one must return to the original derivation of cosmological perturbation theory. At linear order Si(k) = D <?i(k) where 
gi(k) is the normalization and D — a is the growth factor (assuming an EdS Universe as usual), and I have dropped 
the decaying mode solution as usual (I will discuss this further below). At 2nd order 

S 2 (k) = D g 2 (k) + D 2 J -^_ 5l (q) 5l (k~q)4 2) (q,k-q) , (22) 

where the D g 2 (k) term is set to zero in the standard calculation (see [68| for explicit forms of and Jg be- 
low). Again, this equation is missing terms which grow like D^ 1 / 2 or slower (naively one would expect them to be 
unimportant). 

^3 is needed to compute the power spectrum, but computing the bispectrum at this point is informative: 
(5(k 1 )5(k 2 )5(k 3 )> = D 3 (g 1 (k 1 )g 1 (k 2 )g 1 (k 3 )) + D 3 (g 2 (k 1 )g 1 (k 2 )g 1 (k 3 )} 

/73 
- — -3 Jg 2) (q, ki - q) (pi(q)gi(ki - q)gi(k 2 )g 1 (k 3 )) + cyclic permutations of ki .(23) 
(2tt) 

The last term here is the usual perturbative result. The first term is normally zero, for Gaussian initial conditions. 
The middle term is normally zero because g 2 (k) is chosen to be zero. One can look at Eq. (f2"3"]l two ways: The easy 
way is to say that there is no breakdown in perturbation theory, where a higher order term becomes larger than a 
lower order term. The last term is the leading term, i.e., the bispectrum is fundamentally O (Sf), and there is no 
need for the first two terms to enter this discussion at all. The second approach is to argue that all possible terms 
should be included in the calculation - the fact that the first term happens to be zero initially is a peculiarity of the 
initial conditions, not fundamental to the calculation. In this approach, the new middle term, since it is arbitrary, 
could be set in a way that would cancel the last term at some chosen D+. The first term becomes a function of D+, 
and requiring that the solution be independent of -D* would produce an RG equation for (<7i(ki)gi(k2).gi(k3)) (-D*), 
just as above in the simple example. The second approach is clearly more thorough, however, in this paper I will take 
the easy route and declare the usual bispectrum term to be the leading order, with no need for renormalization, i.e., 
I will set g 2 (k) = as usual (we will see that this term is not needed to renormalize the power spectrum). While my 
approach may not give the best possible result, it is not fundamentally incorrect, at least not any more than standard 
perturbation theory is fundamentally incorrect. Renormalization is a tool, not an imperative, and renormalizing part 
of the calculation but not all of it should be better than nothing. 

Moving on to £3, and again including the homogeneous solution, gives 

S 3 (k) = Dg 3 (k) + D 3 f f qi 3 f q t 2 3 gi(qi)gi(q 2 )g 1 (k-qi -q 2 )4 3) (q l ,q 2) k-q 1 -q 2 ) , (24) 

J (27T) (27T) 

where terms that grow like D 1 / 2 or slower have as usual been dropped. The power spectrum is now 

(*(k!)a(k a )) = D 2 (g 1 (k 1 )g 1 (k 2 )) + D 2 ( ff3 (ki) ffl (k 2 ) + g 1 (k 1 )g 3 (k 2 )) + D i (2tt) 3 S d (k x + k 2 )[^]( fcl ) (25) 
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where [<7i](fc) = P13 (A;) + P22(k) is, as [P|](fc) above, just a quick way of writing the usual O (5f) perturbative 
correction terms (note that if I had promoted the leading order bispectrum to O(Sf) there would be more new terms 
in the power spectrum). Following the simple example above, we know just what to do here: eliminate the O (6fj 
term at D+ by a specific choice of 33 (k), i.e., 

( 53 (ki).gi(k 2 ) + ^(ki^M = -Dl(2n) 3 S D (k 1 + k 2 )[.^Pi) , (26) 

which finally leads to 

P(k) = D 2 P n (k) + D 2 (D 2 - Dl)[P 13 (k) + P 22 (k)} , (27) 

where Pn(fc) must be a function of D±. After the same renamings, this is exactly the starting point for the original 
RG calculation above, Eq. (fT3")) . 

While the approach described here apparently takes care of the relative divergence of the corrections to the power 
spectrum of density fluctuations, there is a hidden approximation that guarantees some imperfection. In the usual 
calculation the initial conditions are specified by simply specifying the density-density power spectrum P ss (k), because 
for growing mode initial conditions the power spectrum of velocity divergence fluctuations 9 = -H^V • v, P 9e (k), 
and the cross-correlation between 5 and 9, P 5e , are equal to P ss (k). However, the non- linear corrections to these 
three spectra are not equal, i.e., decaying modes are created. To simultaneously cancel the corrections to all three 
spectra, decaying modes would need to be included in the solutions for 6, e.g., 5% = gi(k)D + di(k)Z?~ 3 / 2 . They do 
not decay instantly, so the non-linear terms they produce can not be safely dropped. It is hard to estimate how large 
the net effect of these changes would be without actually doing the full calculation, although [3!| showed that the 
decaying mode contribution is a modest, although not negligible, fraction of the non-linear evolution (in a somewhat 
different context). Extension of the calculation in this paper to properly account for decaying modes should not be 
too difficult. One simply needs to write down solutions for P 5S (k), P ee (k), and P se (k), including the possibility of 
decaying modes in the initial conditions. This will produce terms similar in form to Pis(k) and Pzvik), with various 
changes in numerical factors. It will then be possible to simultaneously cancel all of the corrections to P ss (k), P 9e (k), 
and P de (k), and obtain RG equations for each of them. If the basic principle behind this method is useful, doing the 
calculation more carefully should result in even better accuracy. 

While the potential usefulness of this RG method appears to be clear, and versions of it have been successful in 
many previous applications [H Q H| US HH [H [H, \M, HE [H [H [H, [Si, [H , the exact conditions in which it 
will succeed or fail are not completely understood. However, the presence of an attractive fixed point (or manifold) 
is clearly a positive feature [49|, [B(j. Eq. (fT4|) has such a fixed point at n — —1.4, which appears to be attractive 
(I haven't proved this but it works for every power spectrum I've tried). In the asymptotic (long time) limit the 
arbitrarily many degrees of freedom in a general power spectrum are reduced to a single amplitude of this power law. 
The system is naturally driven towards a regime where the perturbation theory approximation works best, i.e., the 
2nd order correction to the linear theory power spectrum vanishes. Furthermore, the correction to the bispectrum 
also vanishes for n = —1.4, so we are driven to a place where the Gaussianity of the field is well-preserved. On the 
other hand, in the simple example of Eq. (|15p we obtained an exact solution despite the fact that there is no fixed 
point - in fact, the relative size of the non- linear term diverges with time. Ultimately, the best way to see whether 
this rearrangement of the perturbation series is useful or not is probably to simply compute higher order terms to see 
if the result converges better than standard PT or not. 



III. RESULTS 



In Figs. [HQ] I show the results obtained by solving Eq. (fT4")l . As discussed above, I am assuming an Einstein-de 
Sitter Universe for the background time evolution, to avoid even small time dependence of the [P|] (k) term. I 
compare to the simulation-calibrated fitting formulas of Q , including both their implementation of the formula from 
Peacock and Dodds (hereafter PD96) [7l[ (note that they did not redetermine the parameters using their expanded 
set of simulations) and their newer HALOFIT formula. These formulas were not calibrated for models with baryonic 
acoustic oscillations, so I compare using their standard transfer function from [72| . I choose parameters to produce a 
power spectrum similar to recent best fits to observations as = 0.85, n = 0.96, and T = 0.15. I reiterate that this 
r is chosen to give the closest possible match to the best fit power spectrum from [![, i.e., it has nothing to do with 
the Einstein-de Sitter background evolution, and is less than £l m h in the concordance model because of broad-band 
baryonic power suppression which is not otherwise included in the transfer function of [72j]. 

FigureQJa) shows essentially perfect (better than percent level) agreement at z = between RGPT and the Peacock 
and Dodds fitting formula [71| for k < 0.3 h Mpc -1 . This does not mean that RGPT is perfect - it is likely that 
the agreement with PD96 is partially a coincidence, because we have no reason to expect their formula itself to be 
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FIG. 1: Thick lines show P(k)/P L (k) at z = for the Peacock and Dodds (PD96) fitting formula [U (black, solid), H's 
HALOFIT formula (blue, long-dashed), the RG-improved PT of this paper (red, short-dashed), and standard PT (green, 
dotted). Thin lines show the other cases divided by the prediction of PD96. (a) and (b) are similar except for the axis scales. 
Note that standard PT was used to calibrate HALOFIT at k < 0.1 ftMpc" 1 . 



accurate to this level. In fact, because of the approximation discussed above related to decaying modes, RGPT must 
inevitably have some error as well. I note here that [H's HALOFIT formula is not suitable for testing perturbation 
theory at low-fc because the fit was actually constrained by PT calculations for k < 0.1 /iMpc -1 (R. E. Smith, private 
communication). The fitting formula is a smooth function, so it is not surprising that it does not agree perfectly with 
standard PT at k < 0.1 /iMpc^ 1 , i.e., if the PT power is generally greater than the simulation power, the region near 
k = 0.1/iMpc -1 will be some average of the two, because a step function is not allowed. It is not surprising that 
this imperfect fit would be missed by [H, because the deviation is < 5%, i.e., smaller than their claimed accuracy. 
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Additionally, the k < 0.1 /iMpc -1 power was not constrained exclusively by PT, because the criterion for using scale 
free simulation results was different (R. E. Smith, private communication). In any case it is clear that RGPT is an 
improvement over standard PT. 

Figure |Tf b) shows the comparison at z — again, with rescaled axes, this time focusing on higher k, including 
deep in the non-linear regime. I note optimistically that, while the prediction is not precise at high fc, it continues 
to be a dramatic improvement over linear theory, i.e., even when the non- linear power is a factor of ~ 10 greater 
than the linear power, the perturbation theory result correctly accounts for most of the difference (admittedly, this 
was also true of standard PT). If similar improvements are obtained by removing approximations or simply with 
each additional order of perturbation theory, one could imagine converging quickly to a precise result. However, it 
is possible, or even likely, that missing higher moments of the velocity distribution will undermine the calculation at 
high k. 

One of the notable features of the results is their convergence at high k to the din P/ dink ~ — 1.4 fixed point of 
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Eq. (fT4]l . where the 2nd order PT term vanishes for a pure power law power spectrum [ijj. At z = this limit is 
reached for k > 0.3 /iMpc -1 , while at earlier times it is only reached at higher k (e.g., k > 2>h Mpc -1 at z = 3). 
Qualitatively, this seems like a positive development. Whatever small-scale structure exists initially, it appears to be 
effectively taken out of play. This may alleviate one common worry about this kind of PT, that one is integrating over 
high-fc modes that are highly non-linear, violating the premise of the perturbation theory. This fixed-point behavior 
provides another reason to hope that higher order calculations can significantly extend the range of accuracy of the 
calculation: it may be that the high k behavior of P(k) can be thought of as a slow rolling with scale, as higher 
order terms become important, of this kind of fixed-point power law. The fixed point will also be modified by the full 
inclusion of decaying modes. 

Figure [2] shows the comparison at z — 1, a typical redshift for planned galaxy redshift surveys aimed at measuring 
the baryonic acoustic oscillation feature [22j. As at z = 0, RGPT agrees well with the fitting formulas in the mildly 
non-linear regime, although now the agreement is better with HALOFIT than PD96 at the high k end of this regime. 
An optimistic interpretation of this figure would be that PD96 is more accurate at low k, because it was not constrained 
by standard PT, but HALOFIT becomes more accurate at higher k, as the authors [3] claim it should be. Dedicated 
numerical simulations and/or higher order PT calculations testing convergence are needed to determine the true 
power. Figure [5Jd again shows that the RGPT prediction traces the fitting formulas reasonably well deep into the 
highly non-linear regime. In this case, the agreement is actually as good as the agreement between the two fitting 
formulas, out to k ~ 5/iMpc -1 . 

Moving on to z = 3, relevant to the Lya forest [zll> [zH [Zll , and again potentially future galaxy redshift surveys [22|, 
we see clear deviation between the PT predictions and the fitting formulas, even in the weakly non-linear regime. This 
is especially puzzling considering that the RG and standard version of PT actually give very similar predictions in the 
weakly non-linear regime. This agreement means that the higher order terms that are being effectively re-summed by 
the RG calculation are not important. Considering the substantial disagreement between the two fitting formulas at 
very low k, I am not ready to conclude that PT is wrong. It may be that the fitting formulas were not well-calibrated 
in this regime. Focused simulations are needed to test this. Deeper in the non-linear regime at z = 3, Fig. [Hb) 
shows that the effect of renormalization becomes substantial, and that RGPT agrees quite well with HALOFIT, while 
both disagree significantly with PD96. The optimistic reading is that RGPT is very successful, but considering my 
willingness to dismiss HALOFIT in the weakly non-linear regime, it is probably best to wait for more simulations 
and/or higher order PT calculations before becoming too pleased with RGPT. 

Finally, we consider z = 6. In the weakly non-linear regime shown in Fig. EJa), RG and standard PT are essentially 
identical, and agree quite well with PD96. HALOFIT differs substantially, so again the results are ambiguous. The 
same can be said for the non-linear regime in Fig. |U(b), where the fitting formulas disagree substantially. 

IV. DISCUSSION AND SPECULATION 

The central result of this paper is Eq. (fTl)) , the RG equation for the renormalized power spectrum, demonstrated 
in Fig. []Ja) to give good results for the power spectrum in the quasi-linear regime at low redshift, where standard 
PT performs poorly. It is clear from Figs. [TJHthat the RG improvement to standard PT is generally helpful, in both 
the weakly and strongly non-linear regimes. Except in the most challenging case of z = and high k, the simulation- 
calibrated fitting formula predictions of Peacock and Dodds [7l[ and HALOFIT Q are sufficiently contradictory that 
it is not completely obvious that they are more reliable than RGPT. 

It is possible that this method will turn out to be a computational curiosity, maybe leading to improved a posteriori 
understanding of what we see in simulations, but unable to achieve sufficient accuracy for practical problems, and 
thus leading to little fundamental change in how we carry out precision cosmology measurements. (It may appear 
that the method is guaranteed to be useful for describing weakly non-linear galaxy clustering, but this could still 
be foiled by bias and redshift-space distortions, which I discuss below.) On the other hand, I hope that this line 
of work is only in its infancy, and will lead to a diversification of computational methods, beyond standard N-body 
simulations. I am distinguishing here between systematic approximation methods and heuristic ones like the halo 
model, which rely heavily on calibration by simulations. One of the implicit goals of this work has been to avoid any 
approximations where the path to improved accuracy is unclear. A good quality for an approximation to have is a 
clear limit in which it can be expected to return to the exact result. Note that standard simulations arc themselves 
exactly the kind of systematic approximation that I would like to see more of. Simulating a finite volume with a finite 
number of particles and limited force resolution are all approximations, but each can be tested for convergence as 
their controlling parameters are changed. 

With the hope of stimulating further work, I now discuss a variety of possible directions that could be pursued, 
although I can't guarantee that all of these are good ideas: 

The first priority should be to include the velocity-velocity and velocity-density power spectra as independent 
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FIG. 2: As Fig. [T] except at z ~ 1. We label this and subsequent z > figures by the redshift at which the amplitude of the 
linear theory power in a flat ACDM model with f2 TO = 0.281 matches the amplitude in the calculation. The actual redshift in 
the EdS model used for the calculation is somewhat lower (e.g., for equal z — linear power, z = 1 in the realistic ACDM 
model corresponds to z = 0.61 in the EdS model). 



functions to be renormalized, as I discussed extensively in §11 C 21 The main requirement for doing this properly is 
a straightforward generalization of the standard PT results to allow for general initial conditions, including decaying 
modes. The n ~ —1.4 fixed-point power spectrum will be modified, hopefully to a more interesting, accurate fixed- 
point shape. 

After that, brute force computation to higher order in perturbation theory should produce more accurate results. 
The fact that the current calculation is a substantial improvement over linear theory at all k (as opposed to diverging 
wildly at some point), well beyond the scale where it is no longer very precise, makes me hopeful that each additional 
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order could extend the effective scale significantly. The possibility of estimating the errors in the calculation by 
comparing results at different orders is an equally strong motivation for going to higher order. [In the same sense, the 
present calculation should show where linear theory breaks down much more accurately than the common practice of 
guessing based on the value of A 2 (fc) = fc 3 Pi(fc)/27r 2 .] 

These calculations give clear motivation and targets for high accuracy simulation work. RGPT, being a first 
principles method containing the possibility of internal error control, may in places lead the comparison rather than 
simply being tested itself. If a higher order of PT is computed, and the results agree between different orders to some 
level of precision, it is reasonable to expect that the result is correct, independent of simulations. (Missing higher 
moments of the particle velocity distribution function may be a loophole in this argument, as I discuss below.) With 
the fixed-point concept in mind, it may be interesting to look more carefully at the evolution of the logarithmic slope 
of the power spectrum with time and scale in simulations. 

Abandoning the first principles philosophy, RGPT results may provide a useful template for a new simulation fitting 



12 



1.2 - 



Q 

PL 



& 1 



0.8 - 



T 1 1 1 — I I I 



z~3 








(a) 



J I I I I I L 



1 



J I I I I I L 



0.01 



0.1 
k [h/Mpc] 



FIG. 3: As Fig. [Q except at z ~ 3 (z = 2.08 in the EdS model). 



formula. Because the errors in RGPT are not too large, and should be slowly varying functions of k and cosmological 
parameters, the new formula could be just a very simple parameterization of corrections to RGPT computed using 
simulations. 

One can obviously apply the techniques of this paper to higher order statistics such as the bispectrum (7f| [z3, , 
trispectrum, etc. [79| ■ As discussed in All C 2\ one can consider promoting the leading order for these statistics to the 
more generally possible O(Sf) and O(Sf) (for the bispectrum and trispectrum, respectively). Bispectrum results can 
be compared to the simulation fitting formula of [80]. If the trispectrum computation is sufficiently accurate, RGPT 
should be useful for computing expected statistical errors on the power spectrum, which is notoriously painful to do 
with simulations [8l| . Similarly, errors on higher order statistics may be computable. 

Mostly for convenience, I have assumed that the only effect of deviations from Einstein-de Sitter background 
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evolution is to change the linear growth factor. This is probably sufficient in the weakly non-linear regime, but if 
RGPT is to achieve much accuracy at higher k, this assumption inevitably must be relaxed. This is clear, even though 
the corrections in standard PT are very small, because the dependence of the power spectrum on the equation of state 
of dark energy, w, at fixed observation-time linear theory power [82f . can not be reproduced under the approximation 
used here. There is simply no way the dependence quantified in [82[ can enter the current calculation. Fortunately, 
even small changes in the right hand side of Eq. (|14p can lead to significant changes in the solution of this non- linear 
equation. It will also be necessary to move beyond the very simple time dependence in Eq. j8j) if one wants to include 
any non-gravitational effects that break this form (e.g., gas pressure). 

To describe direct tracers of gas, like the Lya forest or SZ, gas dynamics could be included in the calculation. 
Even for weak lensing, which traces mass directly, gas physics has been shown to matter at a level problematic 
for high precision experiments (83|. Very rudimentary pressure approximations should be relatively straightforward 
to include 84]. While it may be that very sophisticated implementations of heating and cooling are hopeless, one 
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FIG. 4: As Fig. [T] except at z ~ 6 (z = 4.37 in the EdS model). 



should remember that there is some non-trivial inclusion of non-linearity in this calculation, so this is at least worth 
considering. 

To describe clustering of galaxies, or other unavoidably inexact tracers of density, some concept of bias will need to 
be included in the calculation. I think this may be the area where RGPT is most likely to provide a real fundamental 
improvement over simulations, or at least complement to them. With a large but probably achievable amount of 
computer power, one can simulate the power spectrum of dark matter and its dependence on a relatively small 
number of cosmological parameters more or less perfectly. This will not be possible in the foreseeable future for 
galaxies, so robust, high precision cosmological measurements using them will rely on developing a comprehensive 
understanding of how relatively microscopic galaxy formation details translate into large-scale clustering patterns. 
While the separation of scales probably is not clear enough to make the analogy perfect, this is reminiscent of the 
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classic uses of the renormalization group to understand how complex microscopic theories lead to simple macroscopic 
behavior [ir| . The simplest way to introduce bias is to follow 68, 85| in writing the galaxy density fluctuation, S g , as a 
Taylor series in the mass density, 5 g = bid 1 /il, which leads to an expression for the PT power spectrum of galaxies 
involving a few more terms than the mass power spectrum, but not fundamentally more complicated. This path has 
been followed by [86| . Much more ambitiously, one could write local formation rate equations for a galaxy density field 
and apply perturbation theory to them (something like this was discussed using linear theory in (87[). For example, 
p g = p f{p), with p = —p g to represent depletion of gas, and with both components obeying the usual gravitational 
evolution equations. The mean galaxy density as a function of time would be an output of the calculation as well as 
fluctuations. The goal here would not be to write down a realistic local model and compute exact predictions as much 
as to see what kind of relations between observables can be generically expected. An intriguing possibility is that one 
could show that, regardless of the small-scale galaxy formation details, bias can only take certain universal forms, 
described by a small number of free parameters. Note that the validity of the scale-independent bias assumption has 
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never really been proven even in the perfectly linear regime (see [88| for one attempt). This entire discussion of bias 
of galaxies applies equally well to many other traces of density, e.g., the Lya forest [83 |. 

For many applications, redshift space distortions will need to be included. Like bias, these should be relatively 
straightforward to include in the present RG formulation using the extension of [90|'s linear theory calculation de- 
scribed by [68|. [3(| discusses the imperfections of the [§(| approach, and proposes alternative ideas which may be 
useful. 

It may be interesting, as a computational curiosity, to see if RGPT can give sensible results for n > — 1 power law 
initial conditions. In this case the integrals giving the first correction in standard PT truly diverge, in the sense of 
being infinite rather than just large [III- The integrals can be cut off at wavenumber fc c , but the results then depend 
on fc c , which is arbitrary. One may be able to show that, using the RG as described in this paper, the results will 
converge as fc c is taken to infinity. The small- A* integration of Eq. (TTi|) will need to be done increasingly carefully as 
k c is increased, but hopefully convergence to the n = — 1.4 fixed point will quickly erase all trace of the cutoff. Note 
that k c could affect the normalization of the asymptotic power law, which would spoil this idea (e.g., it may only 
work over some limited range of n). 

One potential limitation in RGPT as presented here is the absence of higher moments of the velocity distribution 
function, e.g., the dispersion in the velocities of particles around their mean velocity at a point. I wrote the Vlasov 
equation ([2]) instead of cutting straight to the hydrodynamic equations to point out this problem, shared by standard 
PT (subsequently discussed by [9l|). Dropping higher moments of the velocity distribution, also known as the single- 
stream approximation, is not really an added approximation in standard Eulerian perturbation theory. This is surely 
known to experts (e.g., [92|). but may not be universally recognized. One can easily include these moments as 
variables in the calculation, for example multiplying the Vlasov equation by PiPj and integrating over p gives (after 
some substitution to remove bulk velocity terms) 

^+2«d* + v k d k a» + a* W + a* k d k v> + g40_±jl^ = Q (2g) 

OT 1 + 

where dk — d/dx k , <7 ij (x,t) = (Sv" 1 oV), g lJ (x, r) = (5v % oV 5v k }, with oV the deviation of a particle's velocity 
from the local mean velocity (note that a <r y term would appear in Eq. [5] if I had not dropped it). Unfortunately, 
the lowest order evolution equations contain only the Hubble drag term, so any initial velocity dispersion will self- 
consistently disappear from the perturbation theory. Furthermore, higher order equations contain no source terms, 
meaning that if the velocity dispersion is initially zero it can never become non-zero. In fact, [93| showed that even 
perturbation theory using the distribution function directly leads to the same result. This is a vexing problem because 
stream crossing is clearly ubiquitous in the real Universe; however, [39| estimated that it only becomes significant at 
k > 1 ft.Mpc~ , so it may not be fatal for intermediate scales. It seems likely that, if these variables really matter, 
there should be some way to recover them using a RG method. The basic idea is that, no matter how small the 
effect is in the bare perturbation theory, if it is dynamically relevant it will grow through a RG equation like Eq. (| 14[) . 
Note that even WIMP cold dark matter has a very small seed velocity dispersion J94|]. The solution is probably to 
simply retain the dispersion variable in the calculation, even though it is decaying, in the same way I have argued 
above that one should retain the decaying combination of S and 9 - because it can be rapidly fed by higher order 
terms. Vorticity, V x v, falls out of the standard PT calculation for the same reason velocity dispersion does 
The standard calculation makes use only of the velocity divergence, V • v. A method that successfully reintroduces 
velocity dispersion may also work for vorticity. 

I started this line of work with the intention of using a significantly different type of renormalization group method, 
based on integrating out (averaging over) small-scale (high-fc) modes while modifying the equations describing evo- 
lution of the large-scale modes in a way that preserves their statistics [47j . This method usually starts with a path 
integral formulation of the problem. [93| presented a path integral approach to large-scale structure using the full 
distribution function, and it should not be too hard to write something similar for the usual moments of the dis- 
tribution function. Velocity dispersion should arise inevitably in this approach, because smoothing a dispersionless 
field obviously produces dispersion. It is less inevitable that this dispersion will actually fundamentally change the 
outcome of the calculation, because the effect of these velocities is already present in the standard PT calculation. 
Among other possibly interesting features, the noise discussed in previous RG approaches to large scale structure 
pfj| . |41| may be generated naturally using this approach (i.e., if high-fc modes are eliminated, the time evolution of 
the remaining modes will no longer be perfectly deterministic). 

Further work could elucidate the connection between this work and the re-summation method of |37|. [38l. Related 
renormalization group procedures, starting from general relativistic perturbation theory, are used by [9a] and [96j to 
compute the back-reaction of inhomogeneities on the global evolution of the Universe. 

Finally, the success of the RGPT calculation in this paper, which leads to what looks like a time evolution equation 
for the power spectrum, suggests that it might be useful to consider the exact time evolution equation for the power 
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spectrum, P(k) oc 2 Re (5^ <5_ky, where 5^ is obtained from the usual evolution equations. The usual argument 

against this would be that (^5^ <5_k^ involves higher order terms like the bispectrum, which then require their own 

time evolution equations, and so on, leading to an infinite hierarchy of equations. Firstly, this could be a useful way to 
look at perturbation theory, especially if we are going to be renormalizing statistics like the power spectrum instead of 
the field itself. Second, one may be able to make progress by simply solving these equations numerically, in the same 
sense that simulations are usually used. One could, for example, truncate the hierarchy by replacing the connected 
spectra at some level by the appropriate functions of lower order spectra given by perturbation theory, and look for 
convergence in the results as higher order terms are added. The advantage of working with statistics of the density 
field rather than the density field itself is enormous. One obtains the desired results directly rather than averaging 
over fluctuations in simulations, where there is usually an unavoidable conflict between the need for large boxes to 
limit finite volume effects like sample variance, and the need for high force and mass resolution to properly compute 
small-scale structure. The fact that we are working with continuous, relatively slowing varying functions, rather than 
wildly fluctuating fields, allows many fewer points to be computed, e.g., we can compute a power spectrum over many 
decades of dynamical range by interpolating between a few tens of logarithmically spaced computation points when 
millions or billions of particles would be needed to cover the same range with a simulation. Tasks that at first glance 
seem arduous, e.g., parameterizing the trispectrum or even higher spectra and setting up their evolution equations, 
should be considered in contrast to the decades of person-power expended on simulations. 
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